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ABSTRACT 

Using a reduced proper motion discriminator, I obtain a sample of 4588 sub- 
dwarfs from the Revised NLTT Catalog of Salim & Gould. The ample statistics 
and low contamination permit much more precise determinations of halo parame- 
ters than has previously been possible. The stellar halo is not moving with respect 
to the Local Standard of Rest (LSR) in either the vertical or radial direction, up 
to uncertainties of 2kms~ 1 . This indicates that either the LSR is on a circular 
orbit or the Sun happens to lie very close to an extremum of the LSR's elliptical 
orbit. Similarly, tentative detections of vertical proper motion of Sgr A* relative 
to the LSR are either incorrect or they reflect real physical motion of the central 
black hole relative to the Galactic potential. The correlation coefficients of the 
halo velocity ellipsoid, which would reflect any possible misalignment between 
its principal axes and the cardinal directions of the Galaxy, vanish to within 2%. 
The halo subdwarf luminosity function peaks at My ~ 10.5 with a full width 
half maximum of about 2.5 mag. 

Subject headings: subdwarfs - stars: luminosity function - stars: kinematics - 
stars: statistics 



1. Introduction 

Samples of nearby halo stars can be analyzed to find the bulk properties of the popula- 
tion: their velocity, spatial, and metallicity distributions, as well as their luminosity function. 
The principal difficulty is obtaining a sample that is large enough to draw statistically sig- 
nificant conclusions while still not being contaminated with disk and thick disk stars, which 
locally outnumber halo stars by a factor ~ 10 3 . 

The most secure method to construct such a sample would be to obtain parallaxes, 
proper motions, and radial velocities for a larger, unbiased sample of stars and then select 
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halo stars based on their space motions. The GAIA satellite would be able to do this, but 
even under the most optimistic projections, its data will not be available for well over a 
decade. 

In the absence of such ideal datasets, most nearby halo samples have been culled from 
catalogs of high proper- motion stars (although there are a few notable exceptions to this rule). 
For example, Dahn et al. (1995) obtained trig parallaxes for about 100 proper-motion selected 
stars, thereby determining their transverse velocities, distances, and absolute magnitudes. 
By rigorous selection on transverse velocity (v± > 260 km s -1 ), they obtained a sample 
that was virtually free of disk and thick-disk contamination. The distances and absolute 
magnitudes then allowed them to measure the luminosity function (LF). Of course, to do so 
they had to correct for the halo stars that were eliminated from their sample (along with the 
unwanted disk stars) by their stringent velocity criterion, and this in turn required a model 
of the halo velocity distribution. 

The best such model up to that date was constructed by Casertano, Ratnatunga & 
Bahcall (1990), who used maximum likelihood to decompose two proper-motion selected 
samples into disk, thick disk, and halo components making use of both photometric and 
proper-motion data. They thereby identified different populations within the data, even 
though individual stars could not generally be unambiguously associated with a specific 
population. In particular, Casertano, Ratnatunga & Bahcall (1990) showed that the like- 
lihood fit was significantly improved by allowing for a third "intermediate" or thick disk 
population rather than just two. The kinematics of the halo when so fit were more ex- 
treme than in the 2-component fit earlier obtained by Bahcall & Casertano (1986) because 
thick-disk contamination was drastically reduced. 

RR Lyrae stars are halo tracers selected on variability rather than proper motion. Esti- 
mates of the RR Lyrae absolute magnitude from statistical parallax automatically yield the 
velocity ellipsoid. While this technique has been applied for almost a century, only in the last 
decade or so has it been realized that the RR Lyrae samples are actually mixtures of thick- 
disk and halo stars. Since statistical parallax uses both radial velocities (whose spectra also 
yield metallicity information) and proper motions, and since it derives distances for all stars, 
full kinematic as well as metallicity information is generally available. In a series of papers, 
Layden (1994, 1995, 1997) both systematized pre-existing data and obtained substantial new 
data, thereby laying the basis for a new statistical parallax solution that clearly separated 
the thick-disk and halo populations using a combination of kinematic and metallicity criteria 
(Layden et al. 1996). Popowski & Gould (1998a,b) and Gould k Popowski (1998) (collec- 
tively PG 3 ) introduced new mathematical methods and on this basis conducted a thorough 
overhaul of the Layden et al. (1996) sample, recalibrating much of the old photographic 
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photometry, incorporating more modern extinctions, identifying suspicious astrometry, and 
developing a new method to incorporate non-RR-Lyrae radial velocities into the analysis. 

Of particular note in the present context, PG 3 were the first to measure five of the 
nine components of the halo velocity ellipsoid: all previous analyses had measured the three 
diagonal components of the velocity dispersion tensor and the component of bulk motion 
in the tangential direction (the asymmetric drift), but had assumed that the off-diagonal 
components as well as the bulk motion in the radial and vertical directions were zero. While 
the PG 3 measurements all turned out to be consistent with zero (in the frame of the Local 
Standard of Rest - LSR), the error bars were tantalizingly close to being able to probe some 
interesting scientific questions. 

The bulk motion of the halo in the radial and vertical directions is more likely to coincide 
with the rest frame of the Galaxy than is the motion of the LSR. The LSR could well be on an 
elliptical orbit, in which case it would be moving towards or away from the Galacitic center 
unless the Sun happened to lie at an extremum of this orbit. Indeed, Blitz & Spergel (1991) 
claimed that the LSR is moving outward at 14kms~ 1 based on radial- velocity measurements 
of gas in the outer galaxy (assumed to be on circular orbits). On the other hand, Metzger 
& Schechter (1994) concluded that the LSR was moving inward at 6.6 ± 1.7 km s" 1 based on 
radial velocities of carbon stars in the outer Galaxy. 

Similarly, if the Milky Way disk is warped, then one would expect the LSR to be 
moving either up or down relative to the Galactic rest frame, unless the Sun happened to 
be at an extremum of the warp. Backer & Sramek (1999) found that Sgr A* is moving 
down at 17 ± 6 km s" 1 relative to the LSR. If the supermassive black hole associated with 
Sgr A* is assumed to be at rest with respect to the Galaxy, then this apparent motion would 
actually be a reflex of the warped motion of the LSR. On the other hand, Reid et al. (1999) 
find that Sgr A* is moving in the opposite direction (although with much larger errors) at 
15 ± llkms" 1 . New more precisements measurements are expected soon (M. Reid, private 
communication 2001). 

For a roughly isotropic ensemble of N s stars, the bulk motion C/j can be measured with 
a precision, 



where Cu is the dispersion in the ith directin and n d is the number of components of the 
velocity measured for each star. For the PG 3 sample, N s ~ 170 and rid = 3, while en ~ 
(160 km s^ 1 ) 2 and c 33 ~ ^Okius" 1 ) 2 . Hence <t(E/i) ~ lSkms" 1 and a(U 3 ) ~ Skins" 1 . The 
PG 3 measurement errors were therefore not quite small enough to probe these interesting 
questions. 
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Of course, measurement of a difference between the bulk-halo and LSR velocities would 
not be unambiguous evidence of LSR motion (Oresme 1377). For example, the angular 
momentum vector of material infalling onto the Milky Way could have radically changed 
between the time of the formation of the halo and disk. The former therefore could in 
principle be rotating in a basically polar orbit (albeit with low Mach number) relative to the 
latter. Hence any sort of relative motion would be intriguing evidence of a non-simple Galaxy 
whose exact origins would have to be sorted out making use of other data and arguments. 

Similarly, the off-diagonal elements of the velocity dispersion tensor (normalized to the 
diagonal elements) c\j could potentially provide evidence of asymmetries of the Galaxy that 
would reflect on its origins. The errors in these quantities are ~ (nrf^/3)" 1 / 2 , or about 8% 
for the RR Lyrae sample. To the best of my knowledge, no one has investigated what might 
cause these quantities to differ from zero, so I do not know whether their consistency with 
zero at the 8% level challenges or confirms any theory. Nevertheless, it seems interesting to 
try to probe the off-diagonal elements at higher precision. 

The status of the LF and local density of the stellar halo are also somewhat controversial. 
Dahn et al. (1995) find that the LF peaks at around M v ~ 12, in qualitative agreement with 
the shape of the LF seen in undisturbed globular clusters (Piotto 1997). However, Bahcall & 
Casertano (1986) and Gould, Flynn & Bahcall (1998) find a roughly flat LF over the interval 
9 < My < 13. While Dahn et al. (1995) and Bahcall & Casertano (1986) both studied local 
stars drawn from proper- motion catalogs, Gould, Flynn & Bahcall (1998) adopted a radically 
different approach: they located stars in Hubble Space Telescope images that were too faint at 
their observed color to be in the disk and so were assigned absolute magnitudes and distances 
based on a halo color-magnitude relation. These stars were generally quite distant (> 3 kpc) 
and therefore perhaps not directly comparable to the local samples. Sommer-Larsen & Zhen 
(1990) had earlier suggested that the stellar halo actually has two components, one roughly 
spheroidal and one highly flattened. (The highly flattened component is not to be confused 
with the thick disk: it is not rotating significantly.) In their model, the two components have 
roughly equal densities at the solar circle. Such a model predicts that the halo density should 
be roughly twice as great in the solar neighborhood as it is at a similar Galactocentric radius 
but a few kpc above the Galactic plane. Indeed, Gould, Flynn & Bahcall (1998) found a 
halo density that was lower than the Dahn et al. (1995) measurements by just this fraction. 
Recently Siegel et al. (2002) have argued on the basis of a sample of 70,000 stars along 
multiple pencil beams that even a 2-component halo model is inadequate to explain their 
star counts. 

The newly released Revised NLTT Catalog (Gould & Salim 2002; Salim & Gould 2002b) 
allows one to obtain a very large and very clean sample of halo stars. The reduced proper- 
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motion (RPM) diagram using the newly obtained V — J colors clearly separates main- 
sequence stars, subdwarfs, and white dwarfs into different tracks (Salim & Gould 2002a) in 
sharp contrast to RPM diagram constructed from the original NLTT (Luyten 1979, 1980). 
Although the first release of this catalog covers only 44% of the sky, it contains more than 
5000 local halo stars, well over an order of magnitude more than have ever been cleanly 
distinguished from disk stars on a star-by-star basis. This sample therefore opens the way to 
a much more detailed study of the local halo population than has previously been possible. 
In its present form, the sample does have some limitations. Since most of its stars lack 
radial velocities and parallaxes, it is not possible to establish the absolute distances or the 
amplitude of the velocity ellipsoid based on the Revised NLTT Catalog alone. Nevertheless, 
the amplitude of the velocity ellipsoid is already known with a precision of about 10% 
from previous studies, and by incorporating this external information one can obtain much 
more precise measurements of the five components of the ellipsoid that are currently poorly 
measured: Ui, U3 and c^. Once the velocity scale is set, the mean distances to the stars 
are also determined, which permits one to measure the LF. The large number of stars in 
the sample therefore offers the hope of probing the bottom of the subdwarf sequence which, 
because of its dimness, is poorly represented in magnitude limited samples. Finally, the 
catalog contains a large number of stars from the Galactic plane to about z ~ Vj_/A*ii m ~ 
350 pc above the plane, where V± ~ 300 km s" 1 is the typical transverse speed seen toward 
the Galactic poles and fj,n m = 180 mas yr^ 1 is the proper-motion limit of the catalog. While 
this distance is short compared to the several kpc's hypothesized as the height of the flattened 
halo component, the large number of stars in the sample may yield a statistically significant 
statement about the presence of a density gradient on these larger scales. 

Maximum-likelihood (ML) analysis is absolutely critical for extracting halo parameters 
from this catalog. For example, since the mean tangential velocity of stars seen toward 
the Galactic poles is ~ 200 km s -1 , one might naively expect that the stars selected toward 
the poles would have, on average, this velocity. However, given the fact that the sample 
is proper-motion limited, for all but the dimmest absolute magnitudes the number of stars 
seen with velocities that are la higher than average (300 km s" 1 ) is (300/100) 3 = 27 times 
higher than the number with velocities that are la lower (100 km s" 1 ). This severe selection 
bias does not directly affect any other parameters. However, it couples through the highly 
uneven (but perfectly known) sky coverage from 2MASS, to indirectly affect essentially all 
other parameters. These effects can only be removed by comparing the predictions of models 
with the observations, as ML does automatically. 

Hence, I begin in § 2 by giving a careful summary of the ML modeling procedure. In § 3, 
I present my results and comment on various aspects of these whose interpretation requires 
caution. Finally, in § 4, I compare my results to previous work and briefly discuss the 
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implications of this comparison. I reserve to the Appendix a somewhat technical discussion 
of the problems in determining the completeness of the Revised NLTT Catalog and the 
impact of this completeness on parameter estimation. 

2. Maximum Likelihood Formulation 

2.1. General Equation 

I use ML to estimate the parameters of the stellar halo. In general, a given dataset is 
described by m observables -2 bs,i, which I collectively denote z^ s . If this m-dimensional space 
of observables is divided into bins of volume YliLi ^ z ohs,i, then the likelihood of detecting 
n k objects in the fcth bin is C k = r" k exp(-r fc )/n fc !, where r k = P k {z^ s ) fL=a Az ohSti , and P k 
is the probability density predicted by a given model. If the bins are now made sufficiently 
small that r fe <C 1, then n k < 1 and hence n k \ = 1. The logarithm of the product of the 
likelihoods from all the bins is therefore, 

ln£ = ^ln£ fc = ^n fc lnr fc - ^V fc . (2) 

k k k 

The last term is just the total number of detections expected in the model, N exp . Since the 
n k in the first term are either or 1, equation (2) can be rewritten, 

Ndct m 

In £ = ^ M^(Cbs) II A **-.*l " iVexp- (3) 
k=l i=l 

where iVdet is the total number of detections. In general, however, the probability density 
is not most naturally written directly as a function of the observables z bs,i) but rather of 
the model coordinates, -2 mo d,i, which are evaluated at the observables. Equation (3) can be 
rewritten in terms of these, 

Ndet m 

\nC = J2 HPk[zSod(^L)]J} - N CXP + N dct ^ A^ obS)J . (4) 
k=i i=i 

where J is the Jacobian of the transformation from the observables to the model coordinates. 
For a given set of observations, the last term is the same for all models and so can be dropped. 

2.2. Jacobian 

In the present case, for each star there are six observables: the angular position on the 
sky (/, 6), the proper motion (/^, /i b ), and the two photometric magnitudes (V, J). There are 
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also six model coordinates: the three spatial coordinates r, the two components of transverse 
velocity vj_, and the absolute magnitude, My. Hence, 



J 



d(r,v ± ,M v ) 



d(l,b,fj,i,fj, b ,V,J) 



r cos b 



d(r, M v 



d(V, J) 



(5) 



To evaluate J, I write r and My as implicit functions of V and J, making use of the 
color-magnitude relation, My = F[(V — J)o]- 



10 



0.2[V-A v (r)-M v ]+1 



Mi 



v 



V -J - 



Ay(r) 



(6) 



Ryj J 

1.38 is the ratio of total to 



Here, Ay(r) is the extinction along the line of sight and Ryj 
selective extinction. Partial differentiation of equation (6) yields the matrix equation, 



F'A' 
Rvj 



dr 

dV 
dM v 

dV 



dr 

dJ 
dM v 

dJ 



f r \ 
\F' -F'J 



whose determinant gives, 
J= lnl0 



cosbr 5 F'Q, Q 



In 10 dA v ( F' 



5 d In r\Ryj 



f— -0 

\Rv.7 ) 



(7) 



(8) 



If the reddening vector were parallel to the subdwarf sequence (F'/R VJ ~ 1), then the 
additional term Q in equation (8) would be negligible. In fact, however, this ratio is roughly 
F 1 1 Ryj ~ 2.6, which means that at typical distances r ~ 300 pc and low Galactic latitudes, 
InQ ~ 0.12, and so it cannot be ignored. Note that the prefactor (0.2 In 10 cos b) is an 
irrelevant constant and can be dropped in practical calculations, so that J — > r 5 F'Q. 



2.3. Model Parameters 

I model the stellar halo distribution as the product of a luminosity function (LF), a 
velocity distribution, and a density profile. In addition, halo stars are assumed to obey a 
linear color-magnitude relation 

My = F[(V - J) ] = a(V - J) + b. (9) 

In principle, one might also assume that this relation has some intrinsic dispersion. However, 
for reasons that I discuss below, I do not include such a parameter in the model. Finally, not 
all halo stars satisfying the selection criteria will be detected. I therefore include in the model 
two parameters (Vbreak, /break) that describe the completeness as a function of apparent V 
magnitude. 

C(V) = 1 (V < 12), 
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C(V) 



(Hreak ~ V) + /brcak(V ~ 12) 
(Vbreak ~ 12) 



(12 < V < Kreak), 



C{V) 




(Hreak < V < 20) 



(10) 



This form is motivated by the fact that NLTT is known to be complete to at least V — 11 
away from the plane, and the Revised Catalog has captured essentially all of these stars 
(Gould & Salim 2002). The completeness falls precipitously at faint magnitudes, V > 18. 
The simplest hypothesis is that it is linear over the intervening magnitudes. 

I model the LF with 13 free parameters, one for each magnitude bin centered at My = 3 
to My = 15. I model the velocity distribution as a Gaussian ellipsoid with 9 parameters: 
three for the bulk motion, three for the diagonal components of the dispersion tensor, 
Qj, and three for the off-diagonal components Cjj (i < j). In practice, I use the normalized 
components of the latter (the correlation coefficients) Cy = c^j '(cuCjj) 1 ^ 2 . It is known that 
the halo velocity distribution is highly non-Gaussian, with a kurtosis that is higher than 
Gaussian in the vertical direction, lower than Gaussian in the radial direction, and roughly 
Gaussian in the direction of rotation (PG 3 ). However, modeling the non-Gaussian character 
of the distribution would be quite complicated, and it is straight forward to show that a 
ML fit of a non-Gaussian distribution to a Gaussian model returns unbiased estimates of 
the first two moments. Hence there is no benefit to modelling the non-Gaussian form of the 
distribution unless one wants to investigate the higher moments of the distribution. Since 
these are not a focus of interest in the current paper, I opt for the simpler Gaussian model. 

I model the halo distribution as falling as a power law with Galactocentric distance R 
and exponentially with distance from the Galactic plane z, i.e. p = Pq(R/ Rq)~ u exp(— k\z\). 
The spatial distribution is therefore described by two parameters, v and k, where k may be 
regarded as the inverse scale height. I adopt R = 8kpc. 

Hence, I begin with 28 free parameters, 13 for the LF, 9 for the velocity ellipsoid, 2 for 
the color- magnitude relation, 2 for the completeness fucntion, and 2 for the density profile. 
However, as I now explain, there is one almost perfect degeneracy among these parameters, 
and therefore one of them must be fixed. If the zero-point b of the color- magnitude relation 
(eq. [9]) is increased by Ab, and all the bulk velocities Ui and dispersions y/cu are reduced by 
10~ 2Afe , then all of the model's predicted proper motions will remain unchanged. The only 
difference will be that a star's absolute magnitude (inferred from its color) will be increased 
by AMy = Ab, and so the inferred density of stars of each resulting magnitude bin will be 
increased by 10 a6Ab . Actually, this scaling remains perfect only in the limit Ay — > 0, but 
since the extinction is quite small, the degeneracy is almost perfect. 
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I therefore fix 



U 2 = -2f6.6kms 



-i 



(11) 



the value measured by Gould & Popowski (1998) for their "kinematically selected" sample 
of halo RR Lyrae stars. Gould & Popowski (1998) also evaluated the velocity ellipsoid for 
a "non-kinematically selected" sample of halo stars, the magnitude of whose Ui component 
is somewhat smaller than given by equation (11). However, as I describe in § 2.5 below, 
the present sample is effectively selected using a combination of kinematic and metallicty 
criteria, just as was true for the Gould & Popowski (1998) "kinematically selected" sample. 
Therefore, the scale of the velocity ellipsoid should be fixed by the "kinmatically selected" 
rather than the Gould & Popowski (1998) "non-kinematically selected" sample. 

Finally, as noted above, I do not include a parameter for the dispersion in the color- 
magnitude relation despite the fact that the halo is known to contain a range of metallicities, 
and therefore a range of absolute magnitudes at fixed color, and hence some dispersion 
a {My). Since I do not include this term, the velocity dispersions found by the ML fit will 
be larger than the true dispersions by 



Hence, in principle, if one knew the cu sufficiently well, one could fix them (or one of them or 
their sum), and fit for a (My)- However, as I will show in § 3, the differences Aca are smaller 
than the present uncertainties in the Cu, and therefore this is not a practical possibility. 
Thus, the fit parameters that I am calling u cu" are actually shorthand labels for cu + Aca. 
This means that the derived parameters will not yield new determinations for the Cu. The 
best one can do is use these measurements to place rough upper limits on cr(My). That is, 
among the nine velocity-ellipsoid parameters Ui, cu, c^, new values will be obtained for only 
five: Ui, U 3 , and Cy. 



Among the six observables (I, b, ///, //&, V, J) only the V magnitude has significant errors. 
The angular coordinates (I, b) are known to < 100 mas, about 6 orders of magnitude smaller 
than the scale on which there are significant density gradients. The proper-motion errors 
are 5.5 mas yr^ 1 or 3% of the NLTT proper-motion threshold. Since the intrinsic dispersion 
in proper motions is of order unity, and since the measurement errors add in quadrature to 
these, these errors are utterly negligible. The V errors are about 0.25 mag (Salim & Gould 
2000). As discussed by Salim & Gould (2002b), these are multiplied by 2.1 when entering 
the RPM and therefore cannot be ignored. Finally, the J errors are typically 0.03 mag. 
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2.4. 



Data Characteristics 
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Even though they enter the RPM with somewhat higher weight (3.1), they are then added 
in quadrature to the V errors, and so are also negligible. 

The errors in the original NLTT proper motions were 20masyr _1 (Salim & Gould 
2002b). These proper motions are not used directly in the evaluation of the likelihood, but 
they do have an indirect effect because they influenced Luyten's determination of which 
stars met his proper-motion threshold of 180 mas yr -1 and so which ultimately entered the 
Revised NLTT Catalog (Gould & Salim 2002; Salim & Gould 2002b). 

2.5. Selection Criteria 

One wishes to select as large a sample of halo stars as possible, while effectively excluding 
stars from other populations. Moreover, one would like to restrict selection to those areas 
of the sky with homogeneous completeness characteristics. To achieve the first goal, I make 
use of the discriminator rj introduced by Salim & Gould (2002b), 

v(V rpm ,V - J, shift) = V rpm - 3.1(V - J) - IA7\ sin6| - 2.73. (13) 

where V rpm = V + 5 log(^) is the RPM. Using the (V, V - J) RPM diagram, Salim & Gould 
(2002b) showed that stars in the range < rj < 5.15 are mostly halo stars. To be conservative, 
I restrict the selection by an additional magnitude on each end and require 1 < 17 < 4.15 (see 
Fig. 3 from Salim & Gould 2002b). Completeness of the original NLTT catalog deteriorates 
significantly in areas south of POSS I, and Salim & Gould (2002b) did not even attempt 
to recover faint NLTT stars in this region because their method cannot be applied there. I 
therefore require 5 > — 32°4. Salim & Gould (2002b) showed that in the Galactic latitude 
interval —0.2 < sin b < 0.3, NLTT completeness of main-sequence stars is severely affected, 
dropping from neighboring zones at higher latitude by a factor ~ 10. While the effect is 
much smaller for subdwarfs (and perhaps negligible for white dwarfs), to be conservative, I 
restrict selection to stars outside this range. Finally, of course, the identifications by Salim 
& Gould (2002b) rely critically on 2MASS, and so have only been carried out for the 47% 
of the sky (57% of the region 5 > — 32.°4) that is covered by the second incremental 2MASS 
release. Hence, the spatial selection function alone is quite complex. This fact, together with 
the large number of observables, implies that great care is required to evaluate the likelihood 
function. 

Note that the discriminator rj (eq. [13]) is effectively a function of both kinematics 
and metallicity. That is, rj increases both with higher transverse velocity and with lower 
metallicity (and so lower luminosity at fixed color). It is for this reason that I said in § 2.3 
that the sample is selected using a combination of kinematic and metallicity criteria. 
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2.6. Likelihood Evaluation 

In order to find the model parameters that maximize the likelihood, one must compare in 
a consistent way the likelihood of observing the data given different sets of model parameters. 
This statement is so obvious that it would appear not worth mentioning. However, achieving 
such consistency is by no means trivial. 

The first term in equation (4) is relatively straight forward to calculate because it 
depends only on differential probability functions that are multiplied together. However, the 
second term, the total number of stars expected to enter the sample for a given model is quite 
complicated in several respects. First, as discussed in § 2.5 the selection criteria themselves 
are complex. Second, to estimate iV eX p requires an integral over nine dimensions: six for the 
model coordinates (r, vj_,My), plus three for the measurement errors (V, /^,nltt, A*6,nltt)- 
Recall that even though I am not making use of the NLTT proper-motion measurements, 
they still enter the likelihood function because they affect the sample selection. 

Integration over more than four dimensions is in general more efficiently carried out 
by Monte Carlo than directly. In the present case, the complexity of the 2MASS coverage 
further reinforces the advantages of Monte Carlo integration. However, such an approach 
poses significant difficulties when comparing likelihood estimates at different locations in 
parameter space: Monte Carlo integration introduces Poisson fluctuations into the evaluation 
of iVexp which are of order the square root of the size of the random sample. While these 
fluctuations can to some extent be suppressed by insisting that all realizations have the same 
sample size, the induced fluctuations remain of the same order. To be certain that these do 
not induce roughness in the ML surface of order unity would require Monte Carlo samples 
of C(10 7 ), which would be computationally prohibiive. 

To counter this problem, rather than directly assembling a separate catalog of fake stars 
for each model, I assemble a single catalog of fake stars and assign the stars different weights 
according to the model. That is, I first choose a baseline model that is reasonably close to the 
final model. For each magnitude bin, I assign an absolute magnitude drawn uniformly over 
this bin, a physical location drawn uniformly from the volume within 1 kpc of the Sun, and 
a transverse velocity drawn randomly from the 2-dimensional Gaussian projected-velocity 
distribution expected in the baseline model. I note the Gaussian probability of each such 
fake star, but do not at this point make use of it. 

Next I determine the observational characteristics of the star. For example, I use the 
color-magnitude relation of the baseline model and the star's distance to obtain the true 
(V — J)o color and Vq magnitude. I draw the error in the observed V magnitude from a 
Gaussian distribution, and redden both the color and magnitude according to a simplified 



- 12 - 



extinction law, 



A v (r,b) = 0.075 | esc 6|[1 - exp(-r sin \b\/h d )\ 



(14) 



where hd = 130 pc is the dust scale height. Similarly, I obtain the true proper motion 
from the distance and transverse velocity and draw the two NLTT proper-motion errors 
from a Gaussian distribution. If the NLTT proper motion exceeds 180 mas yr -1 , and the 
RPM discriminator r\ lies within somewhat expanded bounds (—1 < rj < 6.15) as calculated 
within the baseline model, I accept the star into a master list of fake stars. For each model 
I examine each star on this master list and recalculate (V — J) using the model's color- 
magnitude relation. I accept only stars satisfying the selection criterion (1 < rj < 4.15) as 
calculated within the model. Next, I determine iV ex p by counting all of the fake stars thus 
accepted and assigning each one a weight that depends on the model. The weight is a product 
of factors: (R/ Rq)~ v for proximity to the Galactic center, exp(— k\z\) for distance from the 
Galactic plane, C(V) to take account of completeness, as well as a factor for the luminosity 
function of the star's My bin. In particular, I evaluate the probability of the transverse 
velocity given the model and divide this by the tabulated probability of the baeline model. 
In this way, iV ex p is evaluated stochastically for an ensemble of models, without introducing 
random noise into the relative values obtained for different models. 

For purposes of finding the best-fit model, I use a catalog of fake stars drawn from a 
model that is 100 times denser than the actual stellar halo, so that the stochastic character 
of the fake catalog introduces errors in the parameter estimates that are 10 times smaller 
than the Poisson errors. I use bootstrap to calculate the errors in the parameters, i.e., I 
evaluate the scatter in model fits to 25 realizations of the data formed by drawing randomly 
from the actual data with replacement. Each realization is tested against the same catalog 
of fake data. For this purpose, I use a fake catalog drawn from a model that is 10 times 
denser than the actual stellar halo. 



The best-fit model to the 4588 subdwarfs selected according to the criteria described in 
§ 2.5 has the following characteristics. 



3. 



Results 



3.1. Velocity Ellipsoid Parameters 



For the bulk halo motion relative to the Sun, I find 



E/i = 11.4 ±2.2 km s 



,-i 



U 3 = -5.4 ±2.4 km s 



-i 



(15) 
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in the (outward) radial and (upward) vertical directions. Since the Sun moves relative to 
the LSR at —10.0 ± 0.4 km s -1 and 7.2 ± 0.4 km s" 1 in these directions (Dehnen & Binney 
1998), this implies that the LSR is moving relative to the halo at —1.4 ± 2.2 km s -1 radially 
and —1.8 ± 2.4 km s -1 vertically. That is, both components are consistent with zero. If the 
halo is assumed to be stationary in both directions relative to the Galactic potential, then 
either the deviations of the LSR from a circular orbit must be very small, or the Sun must 
lie close to the extrema of these deviations. On the other hand, if the halo is not stationary, 
then it just happens to have almost exactly the same motion as the LSR, which would be a 
most surprising coincidence. 

The errors in equation (15) include only the statistical errors within the fit and not the 
systematic errors induced by fixing the amplitude of the velocity ellipsoid using the Gould 
& Popowski (1998) value for U 2 = —216.6 km s -1 . However, as I now show, this systematic 
error is relatively small. First, the statistical error of U 2 is 12.5 km s -1 or 6%. This induces a 
systematic error in U\ also of 6%, that is, 0.7kms _1 which is small compared to the statistical 
error. There is a second source of error because, while both the present sample and the Gould 
& Popowski (1998) sample were chosen based on a combination of kinematic and metallicity 
criteria, those criteria are not identical, nor even easily comparable. Hence, the values of U 2 
for the two samples need not be identical. It is difficult to judge the size of this systematic 
error, but it is probably also of order lOkms -1 , i.e., about 5%, and therefore again much 
smaller than the statistical error. The systematic errors scale with Ui and therefore are about 
half as big for the vertical motion as the radial motion. 

The three off-diagonal components to the velocity-dispersion tensor are 

ci 2 = 0.024 ±0.014, £i3 = 0.005 ±0.023, c 23 = -0.004 ± 0.026. (16) 

That is, all three are consistent with zero at about the 2% level. (Because there are five 
velocity-ellipsoid parameters being fit, the 1.7 a "detection" of c 12 cannot be regarded as 
even marginally significant.) 

Finally, the three diagonal components are 

^c u + Acjj = (162.4 ± 1.4, 105.8 ± 1.7, 89.4 ± 1.9) kms^ 1 , (17) 

where Acu is defined by equation (12). This can be compared with values of Ca found by 
Gould & Popowski (1998) for kinematically selected RR Lyrae stars. 

v^~= (171 ±10, 99 ±8, 90±7)kms" 1 , (RR Lyraes), (18) 

If the errors in equation (18) were sufficient small, it would be possible to determine the Acu 
in equation (17) and so characterize a(M v ) (the scatter in My at fixed V — J color). See 
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equation (12). However, given the errors, it is immediately clear that a (My) is consistent 
with zero. To find out what upper bound can be put on a(M v ), I first note that because 
C22/U2 ~ 5, most of the potential information comes from the tangential component. At 
the 1 a level, (AC22) 1 / 2 < 55kms _1 . Hence, from equation (12), a(My) < 0.5. This is 
not a very interesting 1 a limit. Moreover, I have not yet incorporated the statistical or 
systematic uncertainties in the amplitude of the velocity ellipsoid as discussed following 
equation (15). Hence, equations (17) and (18) present a reasonably consistent picture, but 
do not significantly constrain a{My). I do note, however, that the comparison of these two 
equations shows that the choice for normalizing the velocity ellipsoid, U 2 = — 216. 6 km s -1 , 
cannot be off by more than about 15%. If it were, then the fit values for (cu + Acu) 1 ^ 2 would 
also change by 15%, and these would then be inconsistent at high significance with the en 
and C33 as measured for RR Lyrae stars. 

3.2. Halo Profile Parameters 

The halo density is not expected to vary much over the small (~ 300 pc) volume that 
is being probed. As discussed in § 2.3, I therefore model the density profile simply as, 
p = p (R/ R )~ v exp(— k\z\). I find, 

i/ = 3.1±1.0, n = 0.022 ± 0.057kpc~ 1 . (19) 

The estimate of v is consistent with many previous determinations which, because they 
are measured over longer baselines, have much small errors. For example, Gould, Flynn & 
Bahcall (1998) find v = 2.96 ±0.27. The n measurement is quite interesting despite the fact 
that (or rather precisely because) it is consistent with zero. At the 2 a level, this constrains 
the scaleheight to be k,^ 1 > 7kpc. If the local halo were composed of two components, 
one highly flattened and one roughly round, then one would expect the density to fall off 
locally over distances that are short compared 7 kpc. Hence, this result should help constrain 
2-component halo models. 

3.3. Color-Magnitude Relation 

I fit for a color-magnitude relation of the form My = a(V — J) + b (eq. [9]) and find 

a = 3.59, b = 0.69. (20) 

The formal uncertainties on these parameters are very small, of order 0.01 mag. However, 
recall from the discussion above equation (11) that b is completely degenerate with the 
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amplitude of the velocity ellipsoid, which was fixed for purposes of the fit but which actually 
has a statistical uncertainty of 6% (and a comparable systematic error). Hence, the true 
error in b is about 0.2 mag. The total error in a is probably not much larger than the formal 
error. 



3.4. Luminosity Function 

The LF is parameterized by 13 separate 1-mag bins, with centers from V — 3 to V — 15. 
I find (as did Gould, Bahcall & Flynn 1997 when they studied the disk LF) that ML estimates 
of the LF tend to magnify Poisson fluctuations according to the following mechanism. First 
suppose that the true LF has a dip at a certain bin. Observational errors will scatter stars 
from the two neighboring bins into this bin, thus tending to wash out the dip. Hence, a 
ML fit, when confronted by a dip in the observed distribution will tend to accentuate it as 
it reconstructs the underlying (true) LF. Now suppose that the true LF is flat over three 
adjacent bins but because of Poisson fluctuations the central bin is depressed. ML will also 
accentuate this dip in an attempt to reconstruct the "true" LF. Hence, particularly for bins 
with low total counts, ML can introduce structure that is not really present. I handle this 
potential problem by imposing a "roughness" penalty AC = 16(difference/sum) 2 , where 
"sum" and "difference" refer to the sum and difference of the LF in each pair of neighboring 
LF bins. This is equivalent to imposing a A% 2 = 1 penalty when neighboring bins differ by 
35%. Thus, if the data really demand a steep gradient, this penalty will permit one, but it 
will squash spurious gradients. 

Figure 1 shows the resulting LF. This LF is significantly correlated with the completeness 
parameters (see eq. [10]), which are derived simultaneously, 

/break = 43 ± 6%, VWcak = 18.27 ± 0.04. (21) 

In the Appendix, I consider arguments that might lead one to suspect that this estimate of 
/break could be substantially too low. I find that these arguments are not compelling and 
therefore adopt the LF calculated using equation (21) as the best estimate. Nevertheless, in 
order to gain a sense of the possible role of such a systematic effect, I also show in Figure 
1 the LF under the assumption that /break = 65%, the highest value that I consider to be 
plausible. 

Figure 2 compares the derived LF (with /break = 43%) to those of several previous 
measurements of the halo LF: those of Bahcall & Casertano (1986), Dahn et al. (1995), 
and Gould, Flynn & Bahcall (1998), which were all previously compared by Gould, Bah- 
call & Flynn (1997). As explained there, the first two LFs are labelled "BC/CRB" and 
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"DLHG/CRB" to indicate that they have been corrected from the originally published LFs 
to take account of the kinematic selection using the velocity ellipsoid of Casertano, Rat- 
natunga & Bahcall (1990), which is very similar to the ellipsoid of Gould & Popowski (1998), 
and to the one derived for the present sample. The present measurement is in reasonably 
good agreement with the DLHB/CRB determination over the range 9 < My < 14 covered 
by the latter. It disagrees strongly with both the BC/CRB and Gould, Flynn & Bahcall 
(1998) determinations (which are in good agreement with each other). The new measure- 
ment extends over a much wider magnitude range and has substantially smaller error bars 
than any previous determination. 

The LF evaluations at My = 15, and to a lesser extent at My = 14, should be interpreted 
cautiously because they depend sensitively on model assumptions. To understand this point, 
one should consider how ML "thinks" when fitting the LF. To zeroth order, it forms an LF 
in the naive way: by counting the number of stars whose dereddened observed (V — J)o 
color and the color- magnitude relation put them in a corresponding My bin, and dividing 
this number by the total effective volume probed by the survey for stars of that My. There 
are respectively 40 and 18 stars in the final two bins, but only roughly 20 and 6 of these 
are assigned to these bins by the final ML LF fit shown in Figures 1 and 2. What prevents 
ML from assigning much higher densities to the LF at these faint magnitudes? In the next 
brightest bin there are 159 stars. This is about four times larger than in the My = 14 bin 
despite the fact that the LF is roughly the same because the effective volume grows rapidly 
with luminosity at these faint magnitudes. As described above, ML takes this as a zeroth- 
order estimate for the number expected in this bin. It then considers how many of these are 
expected to scatter into neighboring bins because the color errors a(V — J) ~ cr(V) = 0.25 
induce errors in My of aa(V) ~ 0.9. That is, roughly 15% of these 159 stars are expected to 
scatter into the My = 14 bin and a few percent into the My = 15 bin. It is by accounting for 
this scatter, as well as scatter form brighter bins, that ML achieves its final estimate. This 
estimate therefore depends quite senstively on the adopted value of cr(V), which is described 
very simply in the model, but could in principle actually be a function of V or of other 
variables. 

In addition, because of the very small number of detected stars in these final bins, there 
is a potential problem of contamination from non-halo stars. Contamination is not generally 
a problem because, as I argued in § 2.5, the discriminator r\ is limited to regions well away 
from main-sequence stars and white dwarfs. However, the density of the very dim stars on 
the RPM diagram is extremely low (see Fig. 3 from Salim & Gould 2002b), so even the low 
residual level of contamination could play a role. 

To establish the LF at My > 14 more securely, and since the total number of stars 
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in the last two bins is very small (58), the simplest approach would be to obtain V^-band 
photometry for all of them. If the ML result is correct, the majority of these will be found 
to have scattered in from brighter bins. Metallicities and radial velocities from spectra of 
the truly red stars could then resolve issues of contamination. The LF of the final two bins 
would then rest on much firmer ground. 

4. Discussion 

To high precision (roughly 2kms _1 ) the LSR is not moving with respect to the halo 
in either the radial or vertical directions. If the halo itself has no radial motion, the first 
result sharply contradicts the conclusion of Blitz & Spergel (1991) based on gas motions 
that the LSR is moving outward at 14kms _1 . On the other hand, it is reasonably consistent 
with the radial-motion estimate of Metzger & Schechter (1994) based on carbon stars. More 
specifically, I find that the halo is moving at 11.4 ± 2.2 km s -1 relative to the Sun, and they 
find that the outer-Galaxy carbon stars are moving at 15.6 ± 1.7 km s -1 . 

I find that all three off-diagonal components of the velocity dispersion tensor are small, 
within ~ 2% of zero. The only previous measurements of these quantities (PG 3 ) were 
consistent with zero, but with errors that were about 4 times larger. To date, I am not 
aware of any effort to predict the off-diagonal terms from theory. 

The measurement presented here of the LF confirms the basic peaked shape found by 
Dahn et al. (1995), but with about 40 times more stars and therefore covering a magnitude 
interval that is roughly twice as large. It is inconsistent with the flat LF found by Bahcall & 
Casertano (1986), and Gould, Flynn & Bahcall (1998) (although in principle, since the latter 
determination was based on stars away from the solar neighborhood, it cannot be rigorously 
ruled out by my measurement). The present measurement is in rough agreement with that 
of Bahcall & Casertano (1986) at brighter magnitudes, My < 9. 

A shortcoming of the present approach is that there is no information about distances 
within the data set, so the scale of the velocity ellipsoid must be set by external information. 
The distance scale could be set by obtaining either radial velocities (RVs) or trig parallaxes 
for a representative (i.e., random) subset of the stars in the sample. The former would yield a 
statistical parallax solution. I stress "representative" because if the subsample is biased, for 
example is weighted toward stars with extreme kinematics and/or low metallicities, then the 
scale of the velocity ellipsoid will be overestimated by statistical parallax because the stars 
with RVs move faster than the those in the sample as a whole. It would be misestimated by 
trig parallax both because the selected stars would be faster and subluminous compared to 
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the sample as a whole. Hence, one must choose a fair sample, and then make use of archival 
data only for stars within that sample. 

From the standpoint of maximizing the precision of the distance-scale measurement 
with the minimum effort, statistical parallax is to be much preferred over trig parallax. 
Even velocity errors of ~ 20kms~ 1 are quite adequate for a statistical parallax measurement 
with the limiting precision o{j])jr\ = 0.65iV _1 / 2 (Popowski & Gould 1998a). Here N is the 
number of stars in the statistical parallax sample and r\ is the distance- scale parameter. 
Even assuming perfect parallaxes, the limit for the trig parallax technique is o~(r))/r] = 
0.2 In 10cr(My)A r_1 / 2 . Given the large number of halo stars in the Revised NLTT Catalog 
and the relative ease of making RV compared to trig parallax measurements, the modest 
per-star advantage of trig parallax will be overwhelmed by the mass-production techniques 
available for RVs. However, good trig parallaxes would provide information on the luminosity 
of individual stars, which cannot be obtained from statistical parallax techniques. Thus, the 
two approaches are complementary. 

I thank Samir Salim for invaluable discussions about both the content and presentation 
of this paper. Work supported by JPL contract 1226901. 

A. Completeness 

The principal source of systematic errors in this analysis is incompleteness (or rather, 
possible misestimation of the completeness) of the Revised NLTT Catalog. As shown in 
Figure 1, such misestimation can significantly affect the determination of the LF. However, 

I find that it does not affect the estimate of the velocity ellipsoid. 

In principle, catalog completeness could be a function of all six observables, i.e., the 
V and J magnitudes, the proper motion, and the position on the sky. While the relative 
completeness of the Revised NLTT compared to the original NLTT is very well understood 
(Salim & Gould 2002b), there is substantially less information available about the absolute 
completeness of the underlying NLTT. 

The completeness of NLTT can be tested either externally or internally. An external 
check requires an independent search for high proper-motion stars, either over the whole sky 
or some fraction of it. By comparing to the Hipparcos (ESA 1997) and Tycho-2 (H0g et al. 
2000) catalogs, Gould & Salim (2002) concluded that NLTT is nearly 100% complete for V < 

II and for Galactic latitudes |6| > 15°. Even near to the plane, completeness is close to 100% 
for n > 400 masyr -1 . Monet et al. (2000) searched for high proper-motion (/i > 400 masyr -1 ) 
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stars to faint magnitudes toward 1378 deg 2 , and found 241 stars, only 17 of which (their 
Tables 2 and 3) they could not match to NLTT. In fact, two of these 17 are actually NLTT 
stars, namely 58785 and 52890, which correspond to entries 1 and 8 from Table 3. Thus, 
over these surveyed areas, NLTT is 94% ± 2% complete for fi > 400 mas yr -1 . However, the 
Luyten Half-Second (LHS) catalog (Luyten 1979), a subset of NLTT which actually extends 
somewhat below 500 mas yr -1 , is almost certaintly substantially more complete than is the 
NLTT at lower proper motions, so this measurement cannot be regarded as representative 
of NLTT as a whole. 

Unfortunately, there are no sytematic studies comparing NLTT detections with an in- 
dependent search for proper motion stars in the range 180masyr _1 < fi < 400masyr _1 and 
at faint magnitudes. In the absence of such external checks, Flynn et al. (2001) conducted 
an internal completeness determination, whose approach is very closely related to the com- 
pleteness measurement carried out here by ML. For a complete sample drawn from a stellar 
population that is uniformly distributed in space, the number of stars in the phase-space 
volume [yUi,/i2] x [Vi, V2], should be (up to Poisson statistics) exactly 8 times the number in 
the volume [2/j.i, 2/^] x [Vi — 5 log 2, Vi — 5 log 2]. This is because the former physical volume 
is 8 times larger while the physical velocities being probed are exactly the same. Flynn et 
al. (2001) made a series of such comparisons at 0.5-magnitude intervals, and by multiplying 
these together, found that the completeness at -Rnltt = 18.5 is 65% that at Rnltt = 13 
(which latter they assumed to be 100%). 

Monet et al. (2000) pointed out that any effect that reduced detections of more distant, 
slower moving stars (relative to faster nearby ones) could masquerade as incompleteness of 
fainter (relative to brighter) stars under this test. In particular, they argued that more distant 
stars would, being on average farther from the plane, have reduced density. Even a small 
difference of 4% in mean density per comparison could lead to the "observed" /break = 65% 
completeness at RNLTT,brcak = 18.5, since it would be multiplied together 11 times between 
-Rnltt = 13 and 18.5. That is, (1 - 0.04) 11 ~ 65%. Monet et al. (2000) presented a figure 
showing that when the procedure is carried out on stars at lower Galactic latitude, the effect 
is much reduced. 

This critique is important because the 2-parameter characterization of incompleteness 
that I use here (see § 2.3) in essence embodies the Flynn et al. (2001) method. As noted in 
§ 3.4, I find /break = 43% ± 6% at K>reak = 18.27. This fraction is consistent with the Flynn 
et al. (2001) result at the 1.5 a level when account is taken of the fact that I have assumed 
100% completeness at V = 12 rather than -Rnltt = 13. The magnitude Hreak is slightly 
brighter than the value found by Flynn et al. (2001), but this is to be expected, since the 
Revised NLTT Catalog depends on 2MASS J-band cross identification. 
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However, the Monet et al. (2000) critique cannot account for the result found in equation 
(21) for several reasons. First, one does not expect the density of halo subdwarfs to fall off 
significantly within the volume probed by NLTT. Second, the ML fit actually allows for such 
a fall-off, so that even if the actual density did not satisfy this theoretical expectation, the 
fit would automatically take account of the fall-off. (In fact according to eq. [19] the best fit 
is consistent with uniform density.) Third, Flynn et al. (2001) are unable to reproduce the 
Monet et al. (2000) figure 3 and instead find relatively comparable results when they apply 
their test to high-latitude and low-latitude stars. 

Nonetheless, it remains possible that there is some other effect that reduces the counts 
of distant relative to nearby stars of the same class. One such effect would be incompleteness 
as a function of proper motion, rather than magnitude. Such incompleteness could be 
very pernicious because if the more distant (and hence slower moving) stars were under- 
represented by a mere 4%, then (as outlined above) the effect would be exponentiated and 
would generate a large apparent incompleteness as a function of magnitude. I therefore 
tested this hypothesis by including a proper-motion term for the completeness within the 
ML fit. However, this term did not improve the fit even slightly. One is thus left without 
any plausible explanation for the relative lack of more distant, slower stars, other than 
incompleteness function of magnitude. 

The little direct information we have on the completeness of NLTT at faint magni- 
tudes and relatively low proper motions is reasonably consistent with the ML estimates of 
completeness derived here. Reid (1990) searched for high proper motion stars in a single 
Schmidt field toward the north Galactic pole. He recovered 63 stars from NLTT as well 
as 15 stars that met NLTT selection criteria but were not in NLTT. This appears to cor- 
respond to a mean completeness of 81%. However, three of the 15 have proper motions 
180 mas yr" 1 < fi < 200 mas yr" 1 , i.e., within la (Salim & Gould 2002b) of the NLTT limit. 
Stars that scatter across the selection boundary are already taken into account in the model- 
ing procedure and should not be counted as due to "incompleteness" . Hence, the true mean 
completeness of NLTT in the Reid (1990) survey area is more like 84 ± 5%. If one restricts 
consideration to stars V < Vb rea k (beyond which the Revised NLTT Catalog is more sensitive 
to the incompleteness of 2MASS than NLTT), then /break = 43 ± 6% corresponds to a mean 
completeness of halo stars of 78 ± 3%, in 1 a agreement with the value just derived from the 
Reid (1990) study. 

In brief, there is no strong evidence to challenge the completeness estimate given by the 
ML fit. 
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Fig. 1. — Logarithm of the luminosity function (LF) derived from 4588 subdwarfs from the 
Revised NLTT Catalog (Gould & Salim 2002; Salim k Gould 2002b). The solid curve and 
open symbols with error bars represent the best fit. The dashed curve and crosses represent 
the fit under the assumption that the derived catalog completeness has been seriously un- 
derestimated. The main features of the LF remain the same. The faintest two bins should 
be interpreted cautiously (see § 3.4). 
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Fig. 2. — Comparison of four halo LFs: The original determinations by Dahn et al. (1995) 
(DLHG) and Bahcall & Casertano (1986) (BC) have been rescaled by Gould, Flynn & Bahcall 
(1998) using the velocity ellipsoid of Casertano, Ratnatunga & Bahcall (1990) (CRB). The 
present work confirms the "bump" in the LF found by DLHG at My ~ 11 as well as the 
fall-off toward brighter mags found by Bahcall & Casertano (1986), but with much smaller 
error bars in both cases. 



